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Abstract 

A  delta-pulse  that  emerges  from  the  origin  of  coordinates  in  a  lossy,  viscous,  homogeneous 
medium  is  studied  as  it  propagates  through  the  medium  experiencing  absorption  and  dispersion. 
Elementary  models  for  such  propagation  problems  have  appeared  in  classical  textbooks  (e.g., 
J.D.  Jackson,  Classical  Electrodynamics,  pp.  212-215,  Wiley,  1962).  Generally  speaking,  as 
the  pulse  propagates,  its  amplitude  decreases,  and  its  width  broadens.  A  more  detailed,  exact 
analysis  of  this  problem  emerges  from  the  study  presented  here  in  which  the  pulse  obeys  a 
third-order  partial  differential  equation  of  parabolic  type  in  time  and  one  space  dimension, 
first  obtained  by  Stokes  in  1885.  We  obtain  the  exact  solution  of  the  pertinent  boundary  and 
initial  value  problem  (BIVP)  posed  in  a  rigorous  fashion,  in  which  the  initial  displacement  is 
the  delta-pulse  in  question.  The  resulting  exact  distortion  of  the  pulse  (i.e.,  amplitude-decaying 
and  width-broadening  shape)  emerges  from  a  series  solution  which  we  have  obtained  by  Laplace 
transform  techniques.  The  solution  exhibits  the  expected  smoothing-out  effects  of  dispersion.  A 
portion  of  the  final  expression,  which  contains  a  sum  of  repeated  integrals  of  the  complementary 
error  function,  is  recast  in  power  series  form,  thus  simplifying  the  final  result.  We  also  show 
how  to  obtain  an  approximate  solution  of  the  present  BIVP  by  means  of  the  method  of  steepest 
descents.  This  approximation  agrees  with  the  leading  term  of  the  exact  solution  described 
above.  Further  examination  of  the  dispersion  relation  associated  with  the  governing  PDE  shows 
that  the  kinematic  viscosity  of  the  medium  must  be  related  to  the  pulse  propagation  speed  in 
a  specific  way  to  insure  that  propagation  with  attenuation  actually  occurs.  Under  this  simple 
but  restrictive  condition,  quantitative  details  of  the  distorted  pulse  amplitude  are  illustrated  in 
several  nondimensional  graphs.  These  are  plotted  versus  time  (or  position)  for  discrete  successive 
values  of  the  spatial  position  (or  time) . 


1  Introduction 

The  propagation  of  a  sound  pulse  through  a  viscous  fluid  is  an  important  problem  that  has  been 
often  revisited  since  the  governing  equation  for  the  propagation  was  established  by  Stokes  [1]  in 
1885.  Various,  more  or  less  generalized,  formulations  have  appeared  in  textbooks  [2,  3,  4,  5]  as 
well  as  in  papers  [6,  7,  8].  The  governing  equation  is  a  third-order  partial  differential  equation  of 
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parabolic  type.  In  one  space  dimension,  the  boundary /initial  value  problem  is 


'  V  d^u  d‘^u  1  d‘^u 

dt  dx"^  dx"^  (?  dt‘^  X  >  ,  t  >  , 

'  u{x,0)  =  0,  u{x,0)  =  0,  (^) 

^  ?x(0,t)  =  uo5{t), 

where  u{x,  t)  is  the  fluid  velocity,  c  is  the  sound  speed,  ly  is  the  kinematic  viscosity  coefficient 
responsible  for  the  attenuation,  ii  is  the  acceleration,  and  5{t)  is  the  Dirac  delta  function.  Attempts 
to  solve  this  equation  exactly,  in  analytic  form,  have  been  unsuccessful  [7] .  Of  particular  interest  is 
the  solution  when  the  boundary  condition  is  an  impulsive  excitation  at  the  origin.  If  such  solution 
were  obtained,  then  it  could  be  used,  by  means  of  the  convolution  theorem,  to  find  solutions  for  any 
other  excitation.  This  basic  solution  will  reveal  the  damping  effects  of  viscosity  in  a  quantitative 
way.  Qualitative  analyses  in  the  literature  [9]  have  exhibited  the  anticipated  behavior  in  which  the 
pulse  amplitude  decreases,  and  its  width  broadens,  as  it  propagates. 

We  note  that  parabolic  PDEs  similar  to  this  often  have  precursor  solutions  with  infinite  speed 
of  propagation,  so  that  disturbances  are  felt  immediately  at  all  points  in  the  fluid.  This  property 
allows  for  the  smothing-out  effect  of  dispersion  within  the  context  of  continuum  mechanics.  This 
issue  has  been  covered  in  detail  by  Weymann  in  1967  [10]  and  will  not  be  considered  here.  In 
the  following  sections,  we  derive  the  analytic  solution  obtained  by  operational  methods  as  well 
as  an  approximation  to  the  analytical  solution  based  on  the  method  of  steepest  descents,  which 
will  identify  the  strongest  contribution  within  the  complete  analytic  solution  to  the  value  of  the 
pertinent  contour  integral. 


2  Theory 

A  Laplace  transform  pair  is  here  defined  as 


/(■s)  =  /  f{t)e  dt, 

Jo 

/W  =  2vri/  • 

^l\  L  J  a  —  lOQ 


where 


=  sf{s)  -/(0+). 

We  take  the  Laplace  transform  of  Eq.  1  with  respect  to  time  t  to  obtain 

(1  +  Ts)uxx{x,  s) - 7^u{x,  s)  =  0, 

{t(0,s)  =  uo, 


(2) 

(3) 

(4) 


where  -u(x,s)  is  the  transform  of  u{x,t),  subscripts  denote  partial  differentiation,  r  =  vjc^,  and, 
because  of  the  initial  condition,  we  take  Uxx{x,  0)  =  0.  Since  u  has  the  units  [L^/T],  r  has  the  units 
of  time.  Because  of  the  transform,  uq  now  takes  the  units  of  length.  The  solution  of  this  system  is 

sx 


u{x,  s)  =  uq  exp 


C\/l  +  TS 


(5) 


We  now  take  the  inverse  Laplace  transform  of  this  last  equation  to  obtain 


u{x,t)  =  Uq- — :  /  e 
ZTTl  J (j—ioo 

which  has  an  essential  singularity  at  s  =  —  1/r.  With  the  change  of  variable 

s  =  p-  1/r, 


Eq.  6  becomes 


/  e-*/-  r-+Vr+ioo 

u[x,t)  =uo— — ^  /  e  ^  dp. 

ZTTl  Jcr+l/r— joo 


ZjTTI  Jcr+l/r— joo 

We  now  define  the  nondimensional  variables 


_  X  -  t 

x  =  t= 

CT  T 


in  which  case  Eq.  8  becomes 


^—t  pa-\-l/r-\-ioo  ^ 

u{x,t)  =  Uo—  /  e 

ZTTI  Jcr+l/r— 200 
^—i  /•cr+l/r+200 


dp 


=  ^0  7^  /  Fi{p)F2{p)eP^  dp, 

ZTTI  Jcr+l/r— 200 


where 


Fi(p)=e-"v^,  F2{p)  =  e^l^. 


If  we  expand  F2{p)  in  a  Maclaurin  series,  the  product  of  these  two  functions  can  be  written 


F{p)  =  Fi{p)F2{p) 


-  -2 

X  X 


=  g-5;v^  I  j _ _ L  _ ^ _ ^ _ t _ ^ _ 

\/Ty/p  2\tp  nW'^l'^p'^/"^ 

.  X  e-^v^  e-^v^  x^ 

=  -I - - - 1 - - - h  •  •  •  H - • - 

y/r  2!r  p  p"-/2 


which  can  be  inverted  term-by-term. 

We  note  that  the  first  term  in  Eq.  15  is  Fi{p),  whose  inverse  fi{t)  is  tabulated  [11,  12];  hence, 

X  F  ,  . 

The  inverse  of  the  second  term  in  Eq.  15  is  also  tabulated  [11]: 


Q-Xy/¥p 


T  y/p  }  T 


_ g-2;V(4i), 


Third  and  subsequent  terms  in  Eq.  15  can  be  inverted  using  the  known  formula  [11] 


Q-dy/P 

pnl2 


}=(4*r/-.,.-erfc(^). 


where  (3  =  x-^/r,  and  erfc(2;)  is  a  symbolic  operator  that  denotes  the  mth  successive  integrals  of 
the  complementary  error  function.  We  recall  that 


2  1 

erfc(2;)  =  dt 

Jz 

and  the  recursion  relation  [13] 

roo 

i"'eiic{z)=  /  i'^~^  eiic{t)  dt,  (n  =  0, 1,  2, . . .). 

J  Z 

Hence  the  inverse  of  the  nth  term  in  Eq.  15  is,  in  nondimensional  form, 


nl'i-nl2 


C 


P 


,n/2 


rn! 


2Vi 


(19) 


(20) 


(21) 


The  complete  inversion  is  thus  obtained  by  combining  Eqs.  16  and  17  with  the  sum  of  all  n  >  2 
terms  in  Eq.  21: 


■y-t 


u{x,  t)  =  Uq - 

T 


X 

y/nt 


/I  \  — 7%  / 

( +  1 )  ^  ^(4t)”/2-i  i-2  erfc  ( 

V2t  /  ^^2  ^ 


.2yft 


(22) 


or,  in  nondimensional  form 
u(x,  t)  =  e“* 


-|  —  CXD  Q7T, _ 2 

^(1  +  +  y  -= - zx^y^  erfc  ,  ^ 

^y{n  +  l)  V2^/! 


X 


(23) 


The  power  series  expansion  of  i^eric{z)  is  given  by  [11] 


1 


i-2zY 


2™  Y{k  +  I)r  (i  +  ^)  ’ 


(24) 


which  can  be  substituted  into  Eq.  23  to  obtain 


u{x,  t)  =  e 


-t 


oo  -n 


i/2  (  — l)^(x/\/F)^ 


2v/t  ^y{n  +  iy’  Y{k  +  i)r  ( 


TT 


n—k 

2 


(25) 


The  triple  factorial  growth  in  the  denominator  of  the  double  sum  is  an  indication  of  the  negligible 
value  of  its  contribution  to  the  solution.  In  terms  of  the  original,  dimensional  variables  x  and  t, 
this  result  is 


u{x,t)  =  UqC 


X 


2cyjTTTt 

p-t/r  oo 


4-  -  ) 


t\n/2 


X 


t  ^y{n  +  l)\T)  \cyi)  ^or(A:  +  l)r(^)  Vc\/rt; 


n  oo 

E 


{-i)‘ 


X  \ 


•  (26) 


The  factor  comes  from  the  shift  in  Eq.  7  when  changing  the  variable  from  s  to  p.  Here  it  is 
more  convenient  to  deal  with  and  plot  the  nondimensional  Eq.  25,  since  a  single  set  of  curves  will 
quantitatively  describe  results  for  all  dimensional  values  of  x  and  t. 


We  note  that,  for  n  =  m  +  2,  the  series  in  Eq.  23  can  be  expressed  as 


m=0  ^  ' 

where  z  =  x/(2\/I).  A  convenient  way  to  evaluate  the  series  is  by  means  of  the  recursion  relation 
[11,  13] 

erfc(z)  =  —  erfc(z)  -  -  erfc(z) .  (27) 

2m  m 

All  m  >  1  terms  can  be  expressed  in  terms  of  the  first  two  by  the  above  relation,  where  the  first 
two  are 

erfc(2:)  =  erfc(2;)  =  erfc(2;).  (28) 


3  An  Approximation 

Another  way  to  proceed  with  the  inversion  of  Eq.  6  is  to  use  the  change  of  variable  from  s  to  p 
given  by 

l  +  rs=p^,  (29) 

which  yields 

u{x,t)=uo -  [  exp  —  pdp.  (30) 

TiiT  Jl  it  \  cpj\ 

In  this  integral,  L  is  a  vertical  path  to  the  right  of  the  origin  O,  and  there  are  an  essential  singularity 


at  p  =  0  and  a  branch  point  at  s  =  —  1/r. 

Let  /(p)  denote  the  function  in  the  exponent  of  Eq.  30, 

f{p)  =  i(p2  _  1)  ("i  _  a")  =  ^  (31) 

r  \  cpj  T  \  pj 

where  a  =  xj {ct) ,  from  which  it  follows  that 

T  f\ 

-/(p)  =  2p-«-^,  (32) 

t  p^ 

The  roots  of  f'{p)  are  the  solutions  of 

2p^  —  ap^  —  a  =  0,  (34) 

which  has  no  negative  real  roots.  We  note  that 

/'(l)=2-2a,  f  (a)  =  a{Q^  -  1) ,  /'(oo)  =  oo.  (35) 

There  are  various  cases  of  interest.  Assume  a  >  1,  in  which  case  the  waveform  that  started  at 
t  =  0  has  not  yet  reached  the  spatial  location  x,  since  x  >  ct.  In  that  case,  from  Eq.  35, 

/'(1)<0,  /'(a)>0,  /'(oo)>0,  (36) 


which,  because  of  the  sign  change,  implies  a  root  between  p  =  1  and  p  =  a.  We  denote  that  root 
Po  (Eig.  1).  Similarly,  if  a  <  1,  there  is  also  a  root  between  p  =  1  and  p  =  a,  except  that,  in 


Figure  1: 


Re  p 


The  p-plane  with  the  locations  of  po,  a,  and  the  path  C  of  steepest  descents. 


this  case,  p  =  a  and  p  =  1  would  exchange  places  in  Fig.  1.  In  both  cases,  f”{po)  >  0.  Thus  the 
path  C  of  steepest  descents  is  parallel  to  the  imaginary  axis.  This  path  passes  through  po  and  is 
asymptotic  to  the  imaginary  axis  for  p  Too.  The  portion  of  the  path  that  contributes  most  to 
the  value  of  the  integral  is  that  near  the  saddle  point  po-  The  other  roots  of  f'{p)  =  0  are  complex 
with  negative  real  parts. 

We  now  write  the  complex  number  p  in  terms  of  its  real  and  imaginary  parts, 

p  =  po  +  ia,  (37) 


so  that  Eq.  30  can  be  written 

1 


u{x,t)  =  uo—[  exp  |-[(po  +  «cr)^  -  1]  1- 
TTzr  Jc  It 


ct{po  +  ia)  \ 

We  can  then  expand  and  approximate  the  exponential  to  obtain 

t 


(po  +  i<T)i  da. 


u{x,t)  K.  uq —  /  exp 
vrr  Jc' 


f{Po) - cr  1  + 


a 


Po/i 


da, 


where  a  =  x/{ct).  Thus, 


where 


u{x,  t)  ^  Uq 


Poe 


f(Po)  f+oc 


ITT 


e-^^^da, 


f{Po)  =  -{pI  -  1)  -  — ) 

T  \  PO/ 


(38) 


(39) 


(40) 


(41) 


and  Po  is  the  root  of  f'{p)  =  0  between  p  =  1  and  p  =  a.  Since  the  integral  in  Eq.  40  has  the 
known  value  \/'k j A,  we  obtain  as  the  final  result 


u{x,  t)  ss  Uq 


Poe 


/(po) 


Po 


TTTt  \j  O'  +Pq 


{a  >  1). 


(42) 


The  dependence  on  x  appears  only  through  a,  which  is  required  to  find  po.  This  solution  is  singular 
at  t  =  0,  which  is  the  time  at  which  the  impulse  5{t)  acts. 


Eq.  42  is  also  valid  for  a  <  1.  In  that  case,  C  does  not  cross  the  essential  singularity  at  p 
However,  Eq.  42  is  not  valid  for  a  =  1.  For  this  case,  Eq.  34  implies 

2p3  _  p2  _  I  — 

which  has  one  real  root  po  =  1  and  two  complex  roots  (— 1  ±  It  is  convenient  for  this 

to  use  the  changes  of  variable 

p  =  l  +  v, 

which  will  have  the  effect  of  introducing  a  shift,  and 

t  =  t'  +  -. 
c 

With  these  changes  of  variables,  Eq.  30  becomes 


u{x,t)  =  no -  [  exp 

mr  Jl 

For  n  <<  1,  this  integral  simplifies  to 

u{x,t)  ^  J 

With  the  additional  change  of  variable 


-{v^  +  2v)[t'  +  -- 


c  c(l  +  n)/J 


(1  +  v)  dv. 


r-\-ioo 


exp 


2v 

T 


t'  + 


2xv 


dv. 


Eq.  47  becomes 


1  ycr  f 

“(■’=.«)  -“0  5- 


^/cfp 

^  ~  2  ’ 

1  a/ct  /■+*“ 

exp 


,  /c^  dp 
xp  +  t\  — 

V  rj  ^ 


which  is  in  the  form  of  an  inverse  Laplace  transform  in  x  [14,  12].  Thus, 


where 


k  = 

T 


^-ky/p 


Since  this  transform  is  tabulated  [11],  we  obtain 

u{x,t)  K.  Uo^ 


Atttx 


However,  since  a  =  1  (and  x  =  ct),  it  therefore  follows  that 


u{x,t)  ^  uo-j==  exp 


\ 

4c^Tt J 


=  0. 

(43) 
case 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

(53) 


where  we  have  included  the  factoring  resulting  from  the  “shift.”  This  result  coincides  with 

the  first  term  of  Eq.  26  with  x  =  ct  in  the  amplitude  factor.  Thus,  the  first  two  terms  of  the  earlier 
solution  are  the  strongest  contributors  to  the  inversion  integral  [15]. 


Re  CO 


'  Im  CO 

ky/c^  —  {ykj^Y 

\  '  ' 

^  . 

• 

C02 

Figure  2:  The  co  solutions  of  the  dispersion  relation. 


4  Dispersion  Relation 

If  we  seek  a  solution  of  Eq.  la  in  the  form 

u(x,t)  =  (54) 

where  k  =  co j c,  we  obtain  the  quadratic 

+  {ivk‘^)co  —  k‘^c^  =  0,  (55) 

whose  solutions  are 

coi^2  =  cizk\J  (?  —  {ykj‘i)‘^  —  iuk^  I2.  (56) 

If  c  >  i^/c/2,  the  solutions  are  complex  with  a  negative  imaginary  part,  as  illustrated  in  Fig.  2. 
In  this  case,  there  is  wave  propagation  through  the  medium  with  an  attenuation  factor  e~'^^ 
which  reduces  the  wave  amplitude  with  time.  The  second  possibility,  c  <  z/ /c/2,  is  of  no  interest, 
since  the  co  solutions  are  purely  imaginary,  and  there  is  no  wave  propagation.  This  situation  then 
resembles  ordinary  diffusion. 

Eq.  55  can  also  be  written  in  terms  of  r  =  v  jc?  ^  in  which  case 

co^  =  (/cc)^(l  —  iwr),  (57) 


which  can  be  solved  for  k: 

k  = 


CO 


:\/l  +  i 


ICOT 


Ca/1  +  {uiry 

=  CV201  +  Wr)-^  ^  ^  ^  ~  ^ 

If  we  now  express  this  complex  propagation  constant  in  the  usual  form 

,  CO 

k  = - h  la, 

Cp 

where  Cp  is  the  phase  velocity,  and  a  is  the  attenuation,  we  obtain 


— 


c\/2a/1  +  {cory 


CO 


a  = 


+  (wTp  -  1 


Vl  +  {^tY  +  1  c^/2  \/l  +  (wt)^ 

We  can  also  write  the  attenuation  factor  a  in  terms  of  the  phase  velocity: 

a/1  +  {cotY  -  1 


c„r 


(58) 

(59) 


(60) 


(61) 


a  = 


(62) 


Figure  3:  Contour  plot  of  nondimensional  velocity  vs.  nondimensional  distance  and  time. 


5  Numerical  Results  and  Discussion 

We  have  obtained  both  exact  and  approximate  solutions  for  the  propagation  behavior  of  a  delta 
function  pulse  that  travels  in  a  viscous  and  dispersive  medium  obeying  a  general  boundary  and 
initial  value  problem  governed  by  a  third-order  partial  differential  equation.  The  exact  solution 
was  obtained  by  operational  methods,  where  the  resulting  infinite  series  was  expressed  in  terms  of 
repeated  integrals  of  the  complementary  error  function.  It  was  also  shown  how  to  tranform  this 
series  into  a  rapidly  convergent  power  series. 

The  approximate  solution  was  found  independently  by  the  method  of  stationary  phase.  It  agrees 
exactly  with  the  first  term  of  the  series  solution  found  before. 

We  display  these  results  in  several  plots  involving  the  nondimensional  solution  presented  in 
Eq.  25.  In  Fig.  3  is  shown  a  contour  plot  of  the  first  two  terms  (omitting  the  double  sum) 
of  the  expression  for  the  nondimensional  fluid  particle  velocity  u  (denoted  ubar  in  the  figure) 
vs.  nondimensional  distance  xbar  from  the  origin  and  nondimensional  time  tbar.  In  Figs.  4  and 
5,  we  show  various  cuts  taken  from  Fig.  3.  In  Fig.  4,  nondimensional  velocity  is  plotted  vs. 
nondimensional  time  for  several  values  of  nondimensional  distance.  In  Fig.  5,  nondimensional 
velocity  is  plotted  vs.  nondimensional  distance  for  several  values  of  nondimensional  time.  All 
three  plots  show  that,  as  the  initial  impulse  advances  either  in  time  or  space,  its  peak  amplitude 
decays,  and  its  width  broadens,  as  would  be  expected  qualitatively.  However,  the  formulas  derived 
(along  with  their  corresponding  plots)  represent  a  quantitative  evaluation  of  the  anticipated  effects 
of  viscosity  and  dispersion. 
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Figure  4:  Plot  of  nondimensional  velocity  vs.  nondimensional  time  for  several  values  of  nondimen¬ 
sional  distance. 
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Figure  5:  Plot  of  nondimensional  velocity  vs.  nondimensional  distance  for  several  values  of  nondi¬ 
mensional  time. 


Thus,  we  have  analytically  and  quantitatively  described  the  decay  and  broadening  of  impulses 
propagating  in  viscous  media  as  modeled  by  Stokes’  classical  boundary-initial  value  problem. 
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